This invention relates to a synthetic aperture radar system (SAR) for recovering data representing terrain height and generating an image thereof and, more particularly, to a system with improved capability of unwrapping the interference pattern provided by two co-registered SAR images. In an interferometric synthetic aperture radar system (IFSAR), complex-valued data representing two SAR images of the same terrain from slightly different angles are collected, processed and then registered, wherein the two images are spatially aligned. One component of the complex-valued data will be the phase of the received radar signals. The data representing the two images are combined pixel by pixel to generate an interference pattern from the images representing the pixel-wise phase differences of the registered images. The resulting phase difference data is proportional to the terrain height in modulo form wherein the modulus is the wavelength of the IFSAR system. In other words, the phase difference data of the combined image is the terrain height remainder after subtracting the multiple of the modulus from the amplitude of terrain height.
The unwrapping of the data representing the interference pattern is the conversion of this data into pixel data representing terrain height.
When there are no noise or aliasing effects in the phase data, phase unwrapping can be carried out by a simple and straight forward manner by a simple integration procedure in which the partial derivatives of the phase data are extracted and then integrated along vertical and horizontal lines. In practice, however, noise and aliasing effects are common in IFSAR data and if a simple integration is applied to IFSAR data corrupted by noise, the resulting data and image produced won't even approximate terrain height because the phase inconsistencies caused by the noise are propagated and accumulated in the integration procedure.
Prior to the present invention, phase unwrapping has been carried out by a least squares unwrapping technique. The least squares approach to phase unwrapping obtains an unwrapped solution by minimizing the differences between the discrete partial derivatives of the (wrapped) phase data and the discrete partial derivatives of the unwrapped solution. Given the wrapped phase .PSI..sub.ij defined on a rectangular grid defined by 0.ltoreq.i.ltoreq.M and 0.ltoreq.j.ltoreq.N, an unwrapped solution .phi..sub.ij on the same grid is sought. The partial derivatives of the unwrapped phase data are defined to be the row and column differences EQU .DELTA..sup.x.sub.ij =.PSI..sub.ij -.PSI..sub.i-1,j, .DELTA..sup.y.sub.ij =.PSI..sub.ij -.PSI..sub.i,j-1. (1)
These differences must be computed as wrapped phase differences, in which the value 2.pi. is added or subtracted as necessary to ensure the differences lie in the interval (-.pi.,.pi.]. The differences are defined at the boundaries of the grid by means of the "boundary conditions" EQU .DELTA..sup.x.sub.0j =-.DELTA..sup.x.sub.1j, .DELTA..sup.y.sub.i0 =-.DELTA..sup.y.sub.i1, .DELTA..sup.x.sub.M+1,j =-.DELTA..sup.x.sub.Mj, .DELTA..sup.y.sub.i,N+1 =-.DELTA..sup.y.sub.iN. (2)
The boundary conditions are essential for guaranteeing a unique solution to the partial differential equation which is derived as Eq. (4) below.
The differences between the partial derivatives of the solution .phi..sub.ij and the partial derivatives defined by Eq. (2) must be minimized. Differentiating the sum EQU .SIGMA..sub.i,j (.phi..sub.ij -.phi..sub.i-1,j -.DELTA..sup.x.sub.ij).sup.2 +.SIGMA..sub.i,j (.phi..sub.ij -.phi..sub.i,j-1 -.DELTA..sup.y.sub.ij).sup.2 ( 3)
with respect to .phi..sub.ij and setting the result equal to zero yields the equation EQU (.phi..sub.i+1,j -2.phi..sub.ij +.phi..sub.i-1,j)+(.phi..sub.ij+1 -2.phi..sub.ij +.phi..sub.i,j-1)=.rho..sub.ij, (4)
where the function .rho..sub.ij is defined by EQU .rho..sub.ij =.DELTA..sup.x.sub.i+1,j -.DELTA..sup.x.sub.ij +.DELTA..sup.y.sub.i,j+1 -.DELTA..sup.y.sub.ij. (5)
The values of the solution array .rho..sub.ij at the boundaries are defined by the boundary conditions EQU .phi..sub.-1,j =.phi..sub.1j, .phi..sub.i,-1 =.phi..sub.i1, .phi..sub.M+1,j =.phi..sub.M-1,j, .phi..sub.i,N+1 =.phi..sub.1,N-1. (6)
The classical method for solving Eq. (4) is called Gauss-Seidel relaxation. In this method, one initializes the solution array .rho..sub.ij to zero and then performs the following updates in an iterative fashion: EQU .phi..sub.ij =[.phi..sub.i+1,j +.phi..sub.i-1,j +.phi..sub.i,j+1 .phi..sub.i,j-1)-.rho..sub.ij ]/4. (7)
In the iteration, each value of .phi..sub.ij in the grid is computed in sequence using the current neighboring values of .phi..sub.ij and the known value of .rho..sub.ij computed from the phase data and the process is repeated until convergence occurs. The Gauss-Seidel relaxation is not a practical method of solving Eq. (4) because of its extremely slow convergence.
Due to its particular form, Eq. (4) can be solved directly (i.e., noniteratively) by means of Fourier techniques rather than Gauss-Seidel relaxation. After an appropriate "mirror reflection" operation, the FFT of the array .rho..sub.ij can be computed, the result divided by an appropriate function of cosines, and the inverse FFT of the result computed to yield the unwrapped solution. Alternatively, discrete cosine transforms can be used more efficiently in place of the FFT's since no mirror reflection is required.
A primary disadvantage of least squares phase unwrapping is that it unwraps through the phase inconsistencies rather than around them. This disadvantage can be remedied with the introduction of weighting values. When certain phase values are known to be corrupt due to noise, aliasing or other degradations, the corrupted phase values should be zero-weighted so they do not affect the unwrapping. In practice, there is provided an array of weighting values 0.ltoreq.w.sub.ij .ltoreq.1 that correspond to the given phase data. The least squares problem then becomes a weighted least squares problem, and the quantity to be minimized becomes EQU .SIGMA..sub.i,j w.sup.x.sub.ij (.phi..sub.ij -.phi..sub.i-1,j -.DELTA..sup.x.sub.ij).sup.2 +.SIGMA..sub.i,j w.sup.y.sub.ij (.phi..sub.ij-.phi..sub.i,j-1 -.DELTA..sup.y.sub.ij).sup.2 (8)
where w.sup.x.sub.ij and w.sup.y.sub.ij are defined by EQU w.sup.x.sub.j =min(w.sub.ij.sup.2, w.sub.i-1,j.sup.2) w.sup.y.sub.ij =min(w.sub.ij.sup.2, w.sub.i,j-1.sup.2). (9)
The least squares solution .phi..sub.ij to this problem is defined by the equation EQU w.sup.x.sub.i+1,j (.phi..sub.i+1,j -.phi..sub.ij)-.sub.w.sup.x.sub.ij (.phi..sub.ij -.phi..sub.i-1,j)+w.sup.y.sub.i,j+1 (.phi..sub.i,j-1 -.phi..sub.ij)-w.sup.y.sub.ij (.phi..sub.ij -.phi..sub.i,j-1)=.rho..sub.ij,(10)
where .rho..sub.ij is defined by EQU .rho..sub.ij =w.sup.x.sub.i+1,j .DELTA..sup.x.sub.i+1,j +w.sup.y.sub.i,j+1 .DELTA..sup.y.sub.i,j+1 -w.sup.x.sub.ij .DELTA..sup.x.sub.ij -w.sup.y.sub.ij .DELTA..sup.y.sub.ij. (11)
Equation (10), unlike the unweighted least squares equation (4), cannot be solved directly. It must be solved iteratively.
The classical Gauss-Seidel relaxation method solves this equation by iterating on the following equation: EQU .phi..sub.ij =(w.sup.x.sub.i+1,j .phi..sub.i+1,j +w.sup.x.sub.ij .phi..sub.i-1,j +w.sup.y.sub.i,j+1 .phi..sub.i,j+1 +w.sup.y.sub.ij .phi..sub.i,j-1 -.rho..sub.ij)/v.sub.ij, (12)
where v.sub.ij is defined by EQU v.sub.ij =w.sup.x.sub.i+1j w.sup.x.sub.ij +w.sup.y.sub.i,j+1 +w.sup.y.sub.ij. (13)
As in the unweighted case, the Gauss-Seidel method converges too slowly to be of practical use.
In an article entitled "Robust Two-Dimensional Weighted and Unweighted Phase Unwrapping that Uses Fast Transforms and Iterative Methods" by D. C. Ghiglia and L. A. Romero in the Journal of Optical Society of America, Vol. 11, No. 1, pp. 107-117, there is described a Preconditioned Conjugate Gradient method for solving the weighted least squares phase unwrapping problem. This method is an iterative technique that solves the weighted least squares unwrapping problem defined by equation (10) and speeds up the convergence of the algorithm.